#################################################################
##LAW BREAKING AND LAW BENDING:##################################
## HOW INTERNATIONAL MIGRANTS NEGOTIATE WITH STATE BORDERS#######
## Schwartz, Simon, Hudson, and Johnson## #######################
## October 2020## ###############################################
### Replication Script Manuscript Body###########################
#################################################################

##### PACKAGES #####
library(multiwayvcov)
library(lmtest)
library(dplyr)
library(ggplot2)
library(stargazer)

##### Data: isq #####
## Load CAT data (full)
isq <- read.csv("isq.csv", header=TRUE)

#Note: All analyses include fixed effects for parish and clustered standard errors for enumeration districts. 
#This was done to account for our sampling design, which blocked by parish and clustered by ED

options(scipen=100)
options(digits=3)

### Creating the aspiring migrant subsample ###
potmig <- subset(isq, aspire.abroad.pre>4)



### Table 1 ###

### Full sample models
## Model 1: enter (illegal migration)
enter.ate <- lm(enterlist.pre ~ enterlist.pre.dum + factor(parish), data=isq)
summary(enter.ate)
enter.ate.clust<-cluster.vcov(enter.ate, isq$ed.parish)
enter.ate.mcse <- coeftest(enter.ate, enter.ate.clust)
enter.ate.mcse

## rownames for stargazer
rownames(enter.ate.mcse)[2] <- "Treat"


## Model 3: work (semi-legal migration)
work.ate <- lm(worklist.pre ~ worklist.pre.dum + factor(parish), data=isq)
summary(work.ate)
work.ate.clust<-cluster.vcov(work.ate, isq$ed.parish)
work.ate.mcse <- coeftest(work.ate, work.ate.clust)
work.ate.mcse

## rownames for stargazer
rownames(work.ate.mcse)[2] <- "Treat"

### Aspiring migrant models

# Because the treatment was randomized at the level of the full sample,
# we include our control variables (gender, age, education, income, networks abroad) 

## Model 2: enter (illegal migration)
enter.ate.potmig <- lm(enterlist.pre ~ enterlist.pre.dum + gender + age + education + income.imp + freq.help.net+  factor(parish), data=potmig)
summary(enter.ate.potmig)
enter.ate.potmig.clust<-cluster.vcov(enter.ate.potmig, potmig$ed.parish)
enter.ate.potmig.mcse <- coeftest(enter.ate.potmig, enter.ate.potmig.clust)
enter.ate.potmig.mcse

## rownames for stargazer
rownames(enter.ate.potmig.mcse)[2] <- "Treat"

## Model 4: work (semi-legal migration)
work.ate.potmig <- lm(worklist.pre ~ worklist.pre.dum + gender + age + education + income.imp + freq.help.net+  factor(parish), data=potmig)
summary(work.ate.potmig)
work.ate.potmig.clust<-cluster.vcov(work.ate.potmig, potmig$ed.parish)
work.ate.potmig.mcse <- coeftest(work.ate.potmig, work.ate.potmig.clust)
work.ate.potmig.mcse

## rownames for stargazer
rownames(work.ate.potmig.mcse)[2] <- "Treat"

## Making the table
stargazer(enter.ate.mcse, enter.ate.potmig.mcse, work.ate.mcse, work.ate.potmig.mcse, align=TRUE,
          keep=c("Treat"), column.labels=c("Full Sample", "Aspiring Migrants", "Full Sample", "Aspiring Migrants"), 
          omit.stat=c("LL","ser","f"),  title="Average Treatment Effects",single.row=FALSE, report="csp")


### Figure 1 ###

### Family and enter (illegal migration)
## running the model
enter.family.potmig <- lm(enterlist.pre ~ enterlist.pre.dum + norm.family.rev + enterlist.pre.dum*norm.family.rev
                          + gender + age + education + income.imp + freq.help.net + factor(parish), 
                          data= potmig)

summary(enter.family.potmig)

##generating clustered standard errors
enter.family.potmig.clust<-cluster.vcov(enter.family.potmig, potmig$ed.parish)
enter.family.potmig.mcse <- coeftest(enter.family.potmig, enter.family.potmig.clust)

## Pull out data for marginal effects plots
# coefficients
enter.family.potmig.b2 <- enter.family.potmig.mcse[2,1]
enter.family.potmig.b3 <- enter.family.potmig.mcse[3,1]
enter.family.potmig.b22 <- enter.family.potmig.mcse[22,1] # interaction term
# variances
enter.family.potmig.varb2 <- enter.family.potmig.clust[2,2]
enter.family.potmig.varb3 <- enter.family.potmig.clust[3,3]
enter.family.potmig.varb22 <- enter.family.potmig.clust[22,22]
# covariances
enter.family.potmig.covb2b22 <- enter.family.potmig.clust[2,22]
enter.family.potmig.covb3b22 <- enter.family.potmig.clust[3,22]

list(enter.family.potmig.b2, enter.family.potmig.b3, enter.family.potmig.b22, enter.family.potmig.varb2, 
     enter.family.potmig.varb3, enter.family.potmig.varb22, enter.family.potmig.covb2b22, enter.family.potmig.covb3b22)

##transforming dataframe to graph    
enter.family.potmig.df <- potmig %>% mutate(norm.family.rev = seq(1, 7, length.out=802),
                                            enter.family.potmig.conbx = enter.family.potmig.b2+enter.family.potmig.b22*norm.family.rev,
                                            enter.family.potmig.consx = sqrt(enter.family.potmig.varb2 + enter.family.potmig.varb22*(norm.family.rev^2)+
                                                                               2*enter.family.potmig.covb2b22*norm.family.rev),
                                            enter.family.potmig.ax = 1.96*enter.family.potmig.consx,
                                            enter.family.potmig.upperx = enter.family.potmig.conbx + enter.family.potmig.ax,
                                            enter.family.potmig.lowerx = enter.family.potmig.conbx - enter.family.potmig.ax)


### Family and work (semi-legal migration)
## running the model

work.family.potmig <- lm(worklist.pre ~ worklist.pre.dum + norm.family.rev + worklist.pre.dum*norm.family.rev
                         + gender + age + education + income.imp + freq.help.net + factor(parish), 
                         data= potmig)

summary(work.family.potmig)

##generating clustered standard errors
work.family.potmig.clust<-cluster.vcov(work.family.potmig, potmig$ed.parish)
work.family.potmig.mcse <- coeftest(work.family.potmig, work.family.potmig.clust)

## Pull out data for marginal effects plots
# coefficients
work.family.potmig.b2 <- work.family.potmig.mcse[2,1]
work.family.potmig.b3 <- work.family.potmig.mcse[3,1]
work.family.potmig.b22 <- work.family.potmig.mcse[22,1] # interaction term
# variances
work.family.potmig.varb2 <- work.family.potmig.clust[2,2]
work.family.potmig.varb3 <- work.family.potmig.clust[3,3]
work.family.potmig.varb22 <- work.family.potmig.clust[22,22]
# covariances
work.family.potmig.covb2b22 <- work.family.potmig.clust[2,22]
work.family.potmig.covb3b22 <- work.family.potmig.clust[3,22]

list(work.family.potmig.b2, work.family.potmig.b3, work.family.potmig.b22, work.family.potmig.varb2, 
     work.family.potmig.varb3, work.family.potmig.varb22, work.family.potmig.covb2b22, work.family.potmig.covb3b22)

## transforming dataframe to graph    
work.family.potmig.df <- potmig %>% mutate(norm.family.rev = seq(1, 7, length.out=802),
                                           work.family.potmig.conbx = work.family.potmig.b2+work.family.potmig.b22*norm.family.rev,
                                           work.family.potmig.consx = sqrt(work.family.potmig.varb2 + work.family.potmig.varb22*(norm.family.rev^2)+
                                                                             2*work.family.potmig.covb2b22*norm.family.rev),
                                           work.family.potmig.ax = 1.96*work.family.potmig.consx,
                                           work.family.potmig.upperx = work.family.potmig.conbx + work.family.potmig.ax,
                                           work.family.potmig.lowerx = work.family.potmig.conbx - work.family.potmig.ax)

### interaction plot with confidence intervals based on clustered SEs

## The rug
family.potmig.temp <- potmig$norm.family.rev
a<-runif(802,min=0,max=0.6)
a<-a-0.2
family.potmig.jitter<-family.potmig.temp+a

## The plot
family.potmig.df <- cbind(work.family.potmig.df, enter.family.potmig.df)
colnames(family.potmig.df) <- make.unique(names(family.potmig.df))

family.potmig.plot.combined <- ggplot(data = family.potmig.df) + 
  theme(text = element_text(size=30)) +
  scale_y_continuous("Estimated Support") + 
  scale_x_continuous(limits = c(1,7)) +
  geom_ribbon(aes(x=norm.family.rev, ymin=work.family.potmig.lowerx, ymax = work.family.potmig.upperx), color="grey", alpha=.2) + 
  geom_line(aes(x=norm.family.rev, y=work.family.potmig.conbx, linetype = "Semi-Legal"), size = 1) + 
  geom_ribbon(aes(x=norm.family.rev, ymin=enter.family.potmig.lowerx, ymax = enter.family.potmig.upperx), color="grey", alpha=.2) + 
  geom_line(aes(x=norm.family.rev, y=enter.family.potmig.conbx, linetype = "Illegal"), size = 1)  +  
  geom_hline(yintercept = 0) +
  geom_rug(aes(x=family.potmig.jitter),sides="b") +
  xlab("OK to Break Law to Help Family") +
  ggtitle("Aspiring Migrant Support\n for Unauthorized Strategies") + 
  theme(plot.title = element_text(face="bold", hjust = 0.5))+ 
  theme(legend.title = element_blank())

family.potmig.plot.combined


### Figure 2 ###

### Unfairness and enter (illegal migration)
## running the model

enter.fair.potmig <- lm(enterlist.pre ~ enterlist.pre.dum + fair.jam+ enterlist.pre.dum*fair.jam
                        + gender + age + education + income.imp + freq.help.net + factor(parish), 
                        data= potmig)

summary(enter.fair.potmig)

## generating clustered standard errors
enter.fair.potmig.clust<-cluster.vcov(enter.fair.potmig, potmig$ed.parish)
enter.fair.potmig.mcse <- coeftest(enter.fair.potmig, enter.fair.potmig.clust)

## Pull out data for marginal effects plots
# coefficients
enter.fair.potmig.b2 <- enter.fair.potmig.mcse[2,1]
enter.fair.potmig.b3 <- enter.fair.potmig.mcse[3,1]
enter.fair.potmig.b22 <- enter.fair.potmig.mcse[22,1] # interaction term
# variances
enter.fair.potmig.varb2 <- enter.fair.potmig.clust[2,2]
enter.fair.potmig.varb3 <- enter.fair.potmig.clust[3,3]
enter.fair.potmig.varb22 <- enter.fair.potmig.clust[22,22]
# covariances
enter.fair.potmig.covb2b22 <- enter.fair.potmig.clust[2,22]
enter.fair.potmig.covb3b22 <- enter.fair.potmig.clust[3,22]

list(enter.fair.potmig.b2, enter.fair.potmig.b3, enter.fair.potmig.b22, enter.fair.potmig.varb2, 
     enter.fair.potmig.varb3, enter.fair.potmig.varb22, enter.fair.potmig.covb2b22, enter.fair.potmig.covb3b22)

##transforming dataframe to graph    
enter.fair.potmig.df <- potmig %>% mutate(fair.jam = seq(1, 7, length.out=802),
                                          enter.fair.potmig.conbx = enter.fair.potmig.b2+enter.fair.potmig.b22*fair.jam,
                                          enter.fair.potmig.consx = sqrt(enter.fair.potmig.varb2 + enter.fair.potmig.varb22*(fair.jam ^2)+
                                                                           2*enter.fair.potmig.covb2b22*fair.jam),
                                          enter.fair.potmig.ax = 1.96*enter.fair.potmig.consx,
                                          enter.fair.potmig.upperx = enter.fair.potmig.conbx + enter.fair.potmig.ax,
                                          enter.fair.potmig.lowerx = enter.fair.potmig.conbx - enter.fair.potmig.ax)


### Unfairness and work (semi-legal migration)
## running the model

work.fair.potmig <- lm(worklist.pre ~ worklist.pre.dum + fair.jam + worklist.pre.dum*fair.jam
                       + gender + age + education + income.imp + freq.help.net + factor(parish), 
                       data= potmig)

summary(work.fair.potmig)

## generating clustered standard errors
work.fair.potmig.clust<-cluster.vcov(work.fair.potmig, potmig$ed.parish)
work.fair.potmig.mcse <- coeftest(work.fair.potmig, work.fair.potmig.clust)

## Pull out data for marginal effects plots
# coefficients
work.fair.potmig.b2 <- work.fair.potmig.mcse[2,1]
work.fair.potmig.b3 <- work.fair.potmig.mcse[3,1]
work.fair.potmig.b22 <- work.fair.potmig.mcse[22,1] # interaction term
# variances
work.fair.potmig.varb2 <- work.fair.potmig.clust[2,2]
work.fair.potmig.varb3 <- work.fair.potmig.clust[3,3]
work.fair.potmig.varb22 <- work.fair.potmig.clust[22,22]
# covariances
work.fair.potmig.covb2b22 <- work.fair.potmig.clust[2,22]
work.fair.potmig.covb3b22 <- work.fair.potmig.clust[3,22]

list(work.fair.potmig.b2, work.fair.potmig.b3, work.fair.potmig.b22, work.fair.potmig.varb2, 
     work.fair.potmig.varb3, work.fair.potmig.varb22, work.fair.potmig.covb2b22, work.fair.potmig.covb3b22)

## transforming dataframe to graph    
work.fair.potmig.df <- potmig %>% mutate(fair.jam = seq(1, 7, length.out=802),
                                         work.fair.potmig.conbx = work.fair.potmig.b2+work.fair.potmig.b22*fair.jam ,
                                         work.fair.potmig.consx = sqrt(work.fair.potmig.varb2 + work.fair.potmig.varb22*(fair.jam ^2)+
                                                                         2*work.fair.potmig.covb2b22*fair.jam ),
                                         work.fair.potmig.ax = 1.96*work.fair.potmig.consx,
                                         work.fair.potmig.upperx = work.fair.potmig.conbx + work.fair.potmig.ax,
                                         work.fair.potmig.lowerx = work.fair.potmig.conbx - work.fair.potmig.ax)

### interaction plot with confidence intervals based on clustered SEs

# The rug
fair.potmig.temp <- potmig$fair.jam
a<-runif(802,min=0,max=0.6)
a<-a-0.2
fair.potmig.jitter<-fair.potmig.temp+a

# The plot
fair.potmig.df <- cbind(work.fair.potmig.df, enter.fair.potmig.df)
colnames(fair.potmig.df) <- make.unique(names(fair.potmig.df))

fair.potmig.plot.combined <- ggplot(data = fair.potmig.df) + 
  theme(text = element_text(size=30)) +
  scale_y_continuous("Estimated Support") + 
  scale_x_continuous(limits = c(1,7)) +
  geom_ribbon(aes(x=fair.jam, ymin=work.fair.potmig.lowerx, ymax = work.fair.potmig.upperx), color="grey", alpha=.2) + 
  geom_line(aes(x=fair.jam, y=work.fair.potmig.conbx, linetype = "Semi-Legal"), size = 1) + 
  geom_ribbon(aes(x=fair.jam, ymin=enter.fair.potmig.lowerx, ymax = enter.fair.potmig.upperx), color="grey", alpha=.2) + 
  geom_line(aes(x=fair.jam, y=enter.fair.potmig.conbx, linetype = "Illegal"), size = 1)  +  
  geom_hline(yintercept = 0) +
  geom_rug(aes(x=fair.potmig.jitter),sides="b") +
  xlab("Immigration Authorities are Unfair") +
  ggtitle("Aspiring Migrant Support\n for Unauthorized Strategies") + 
  theme(plot.title = element_text(face="bold", hjust = 0.5))+ 
  theme(legend.title = element_blank())

fair.potmig.plot.combined

### Figure 3 ###

### Perceived risk and enter (illegal migration)
## running the model

enter.risk.potmig <- lm(enterlist.pre ~ enterlist.pre.dum + risk.enter.dest + enterlist.pre.dum*risk.enter.dest 
                        + gender + age + education + income.imp + freq.help.net + factor(parish), 
                        data= potmig)

summary(enter.risk.potmig)

## generating clustered standard errors
enter.risk.potmig.clust<-cluster.vcov(enter.risk.potmig, potmig$ed.parish)
enter.risk.potmig.mcse <- coeftest(enter.risk.potmig, enter.risk.potmig.clust)

## Pull out data for marginal effects plots
# coefficients
enter.risk.potmig.b2 <- enter.risk.potmig.mcse[2,1]
enter.risk.potmig.b3 <- enter.risk.potmig.mcse[3,1]
enter.risk.potmig.b22 <- enter.risk.potmig.mcse[22,1] # interaction term
# variances
enter.risk.potmig.varb2 <- enter.risk.potmig.clust[2,2]
enter.risk.potmig.varb3 <- enter.risk.potmig.clust[3,3]
enter.risk.potmig.varb22 <- enter.risk.potmig.clust[22,22]
# covariances
enter.risk.potmig.covb2b22 <- enter.risk.potmig.clust[2,22]
enter.risk.potmig.covb3b22 <- enter.risk.potmig.clust[3,22]

list(enter.risk.potmig.b2, enter.risk.potmig.b3, enter.risk.potmig.b22, enter.risk.potmig.varb2, 
     enter.risk.potmig.varb3, enter.risk.potmig.varb22, enter.risk.potmig.covb2b22, enter.risk.potmig.covb3b22)

## transforming dataframe to graph    
enter.risk.potmig.df <- potmig %>% mutate(risk.enter.dest = seq(0, 10, length.out=802),
                                          enter.risk.potmig.conbx = enter.risk.potmig.b2+enter.risk.potmig.b22*risk.enter.dest,
                                          enter.risk.potmig.consx = sqrt(enter.risk.potmig.varb2 + enter.risk.potmig.varb22*(risk.enter.dest^2)+
                                                                           2*enter.risk.potmig.covb2b22*risk.enter.dest),
                                          enter.risk.potmig.ax = 1.96*enter.risk.potmig.consx,
                                          enter.risk.potmig.upperx = enter.risk.potmig.conbx + enter.risk.potmig.ax,
                                          enter.risk.potmig.lowerx = enter.risk.potmig.conbx - enter.risk.potmig.ax)


### Perceived risk and work (semi-legal migration)
## running the model

work.risk.potmig <- lm(worklist.pre ~ worklist.pre.dum + risk.work.dest + worklist.pre.dum*risk.work.dest 
                       + gender + age + education + income.imp + freq.help.net + factor(parish), 
                       data= potmig)

summary(work.risk.potmig)

## generating clustered standard errors
work.risk.potmig.clust<-cluster.vcov(work.risk.potmig, potmig$ed.parish)
work.risk.potmig.mcse <- coeftest(work.risk.potmig, work.risk.potmig.clust)

## Pull out data for marginal effects plots
# coefficients
work.risk.potmig.b2 <- work.risk.potmig.mcse[2,1]
work.risk.potmig.b3 <- work.risk.potmig.mcse[3,1]
work.risk.potmig.b22 <- work.risk.potmig.mcse[22,1] # interaction term
# variances
work.risk.potmig.varb2 <- work.risk.potmig.clust[2,2]
work.risk.potmig.varb3 <- work.risk.potmig.clust[3,3]
work.risk.potmig.varb22 <- work.risk.potmig.clust[22,22]
# covariances
work.risk.potmig.covb2b22 <- work.risk.potmig.clust[2,22]
work.risk.potmig.covb3b22 <- work.risk.potmig.clust[3,22]

list(work.risk.potmig.b2, work.risk.potmig.b3, work.risk.potmig.b22, work.risk.potmig.varb2, 
     work.risk.potmig.varb3, work.risk.potmig.varb22, work.risk.potmig.covb2b22, work.risk.potmig.covb3b22)

## transforming dataframe to graph    
work.risk.potmig.df <- potmig %>% mutate(risk.work.dest = seq(0, 10, length.out=802),
                                         work.risk.potmig.conbx = work.risk.potmig.b2+work.risk.potmig.b22*risk.work.dest,
                                         work.risk.potmig.consx = sqrt(work.risk.potmig.varb2 + work.risk.potmig.varb22*(risk.work.dest^2)+
                                                                         2*work.risk.potmig.covb2b22*risk.work.dest),
                                         work.risk.potmig.ax = 1.96*work.risk.potmig.consx,
                                         work.risk.potmig.upperx = work.risk.potmig.conbx + work.risk.potmig.ax,
                                         work.risk.potmig.lowerx = work.risk.potmig.conbx - work.risk.potmig.ax)

### interaction plot with confidence intervals based on clustered SEs

# The rug
enter.risk.potmig.temp <- potmig$risk.enter.dest
work.risk.potmig.temp <- potmig$risk.work.dest
a<-runif(802,min=0,max=0.6)
a<-a-0.2
enter.risk.potmig.jitter<-enter.risk.potmig.temp+a
work.risk.potmig.jitter<-work.risk.potmig.temp+a

# The plot
risk.potmig.df <- cbind(work.risk.potmig.df, enter.risk.potmig.df)
colnames(risk.potmig.df) <- make.unique(names(risk.potmig.df))

risk.potmig.plot.combined <- ggplot(data = risk.potmig.df) + 
  theme(text = element_text(size=30)) +
  scale_y_continuous("Estimated Support") + 
  scale_x_continuous(limits = c(0,10), labels=c(0, 2.5, 5.0, 7.5, 10)) +
  geom_ribbon(aes(x=risk.work.dest, ymin=work.risk.potmig.lowerx, ymax = work.risk.potmig.upperx), color="grey", alpha=.2) + 
  geom_line(aes(x=risk.work.dest, y=work.risk.potmig.conbx, linetype = "Semi-Legal"), size = 1) + 
  geom_ribbon(aes(x=risk.work.dest, ymin=enter.risk.potmig.lowerx, ymax = enter.risk.potmig.upperx), color="grey", alpha=.2) + 
  geom_line(aes(x=risk.work.dest, y=enter.risk.potmig.conbx, linetype = "Illegal"), size = 1)  +  
  geom_hline(yintercept = 0) +
  geom_rug(aes(x=enter.risk.potmig.jitter),sides="t", color="grey") +
  geom_rug(aes(x=work.risk.potmig.jitter),sides="b", color="black") +
  xlab("Perceived Risk") +
  ggtitle("Aspiring Migrant Support\n for Unauthorized Strategies") + 
  theme(plot.title = element_text(face="bold", hjust = 0.5))+ 
  theme(legend.title = element_blank())

risk.potmig.plot.combined

### Figure 4 ###

### Punishment and ener (illegal migration)
## running the model

enter.punishment.potmig <- lm(enterlist.pre ~ enterlist.pre.dum + punishment.enter + enterlist.pre.dum*punishment.enter
                              + gender + age + education + income.imp + freq.help.net + factor(parish), 
                              data= potmig)

summary(enter.punishment.potmig)

## generating clustered standard errors
enter.punishment.potmig.clust<-cluster.vcov(enter.punishment.potmig, potmig$ed.parish)
enter.punishment.potmig.mcse <- coeftest(enter.punishment.potmig, enter.punishment.potmig.clust)

## Pull out data for marginal effects plots
# coefficients
enter.punishment.potmig.b2 <- enter.punishment.potmig.mcse[2,1]
enter.punishment.potmig.b3 <- enter.punishment.potmig.mcse[3,1]
enter.punishment.potmig.b22 <- enter.punishment.potmig.mcse[22,1] # interaction term
# variances
enter.punishment.potmig.varb2 <- enter.punishment.potmig.clust[2,2]
enter.punishment.potmig.varb3 <- enter.punishment.potmig.clust[3,3]
enter.punishment.potmig.varb22 <- enter.punishment.potmig.clust[22,22]
# covariances
enter.punishment.potmig.covb2b22 <- enter.punishment.potmig.clust[2,22]
enter.punishment.potmig.covb3b22 <- enter.punishment.potmig.clust[3,22]

list(enter.punishment.potmig.b2, enter.punishment.potmig.b3, enter.punishment.potmig.b22, enter.punishment.potmig.varb2, 
     enter.punishment.potmig.varb3, enter.punishment.potmig.varb22, enter.punishment.potmig.covb2b22, enter.punishment.potmig.covb3b22)

## transforming dataframe to graph    
enter.punishment.potmig.df <- potmig %>% mutate(punishment.enter = seq(1, 4, length.out=802),
                                                enter.punishment.potmig.conbx = enter.punishment.potmig.b2+enter.punishment.potmig.b22*punishment.enter,
                                                enter.punishment.potmig.consx = sqrt(enter.punishment.potmig.varb2 + enter.punishment.potmig.varb22*(punishment.enter^2)+
                                                                                       2*enter.punishment.potmig.covb2b22*punishment.enter),
                                                enter.punishment.potmig.ax = 1.96*enter.punishment.potmig.consx,
                                                enter.punishment.potmig.upperx = enter.punishment.potmig.conbx + enter.punishment.potmig.ax,
                                                enter.punishment.potmig.lowerx = enter.punishment.potmig.conbx - enter.punishment.potmig.ax)


### Punishment and work (semi-legal migration)
## running the model

work.punishment.potmig <- lm(worklist.pre ~ worklist.pre.dum + punishment.work + worklist.pre.dum*punishment.work
                             + gender + age + education + income.imp + freq.help.net + factor(parish), 
                             data= potmig)

summary(work.punishment.potmig)

## generating clustered standard errors
work.punishment.potmig.clust<-cluster.vcov(work.punishment.potmig, potmig$ed.parish)
work.punishment.potmig.mcse <- coeftest(work.punishment.potmig, work.punishment.potmig.clust)

## Pull out data for marginal effects plots
# coefficients
work.punishment.potmig.b2 <- work.punishment.potmig.mcse[2,1]
work.punishment.potmig.b3 <- work.punishment.potmig.mcse[3,1]
work.punishment.potmig.b22 <- work.punishment.potmig.mcse[22,1] # interaction term
# variances
work.punishment.potmig.varb2 <- work.punishment.potmig.clust[2,2]
work.punishment.potmig.varb3 <- work.punishment.potmig.clust[3,3]
work.punishment.potmig.varb22 <- work.punishment.potmig.clust[22,22]
# covariances
work.punishment.potmig.covb2b22 <- work.punishment.potmig.clust[2,22]
work.punishment.potmig.covb3b22 <- work.punishment.potmig.clust[3,22]

list(work.punishment.potmig.b2, work.punishment.potmig.b3, work.punishment.potmig.b22, work.punishment.potmig.varb2, 
     work.punishment.potmig.varb3, work.punishment.potmig.varb22, work.punishment.potmig.covb2b22, work.punishment.potmig.covb3b22)

## transforming dataframe to graph    
work.punishment.potmig.df <- potmig %>% mutate(punishment.work = seq(1, 4, length.out=802),
                                               work.punishment.potmig.conbx = work.punishment.potmig.b2+work.punishment.potmig.b22*punishment.work,
                                               work.punishment.potmig.consx = sqrt(work.punishment.potmig.varb2 + work.punishment.potmig.varb22*(punishment.work^2)+
                                                                                     2*work.punishment.potmig.covb2b22*punishment.work),
                                               work.punishment.potmig.ax = 1.96*work.punishment.potmig.consx,
                                               work.punishment.potmig.upperx = work.punishment.potmig.conbx + work.punishment.potmig.ax,
                                               work.punishment.potmig.lowerx = work.punishment.potmig.conbx - work.punishment.potmig.ax)


## interaction plot with confidence intervals based on clustered SEs

# The rug
enter.punishment.potmig.temp <- potmig$punishment.enter
work.punishment.potmig.temp <- potmig$punishment.work
a<-runif(802,min=0,max=0.6)
a<-a-0.2
enter.punishment.potmig.jitter<-enter.punishment.potmig.temp+a
work.punishment.potmig.jitter<-work.punishment.potmig.temp+a

# The plot
punishment.potmig.df <- cbind(work.punishment.potmig.df, enter.punishment.potmig.df)
colnames(punishment.potmig.df) <- make.unique(names(punishment.potmig.df))

punishment.potmig.plot.combined <- ggplot(data = punishment.potmig.df) + 
  theme(text = element_text(size=30)) +
  scale_y_continuous("Estimated Support") + 
  scale_x_discrete(limits = c(1,4)) +
  geom_ribbon(aes(x=punishment.work, ymin=work.punishment.potmig.lowerx, ymax = work.punishment.potmig.upperx), alpha=.2) + 
  geom_line(aes(x=punishment.work, y=work.punishment.potmig.conbx, linetype = "Semi-Legal"), size = 1) + 
  geom_ribbon(aes(x=punishment.enter.1, ymin=enter.punishment.potmig.lowerx, ymax = enter.punishment.potmig.upperx), alpha=.2) + 
  geom_line(aes(x=punishment.enter.1, y=enter.punishment.potmig.conbx, linetype = "Illegal"), size = 1)  +  
  geom_hline(yintercept = 0) +
  geom_rug(aes(x=enter.punishment.potmig.jitter),sides="t", color="grey") +
  geom_rug(aes(x=work.punishment.potmig.jitter),sides="b", color="black") +
  xlab("Perceived Severity") +
  ggtitle("Aspiring Migrant Support\n for Unauthorized Strategies") + 
  theme(plot.title = element_text(face="bold", hjust = 0.5))+
  theme(legend.title = element_blank())

punishment.potmig.plot.combined

### Table 2 ###

### Model 1: enter (illegal migration) 

## running the model
enter.potmig <- lm(enterlist.pre ~ enterlist.pre.dum  + fair.jam + norm.family.rev
                   + enterlist.pre.dum*fair.jam +  enterlist.pre.dum*norm.family.rev + risk.enter.dest
                   + risk.enter.dest*enterlist.pre.dum + punishment.enter + punishment.enter*enterlist.pre.dum 
                   + gender + gender*enterlist.pre.dum + age + age*enterlist.pre.dum 
                   + education + education*enterlist.pre.dum + income.imp + income.imp*enterlist.pre.dum 
                   + freq.help.net + freq.help.net*enterlist.pre.dum + factor(parish), data= potmig)
summary(enter.potmig)

## generate clustered standard errors 

enter.potmig.clust<-cluster.vcov(enter.potmig, potmig$ed.parish)
enter.potmig.mcse <- coeftest(enter.potmig, enter.potmig.clust)
enter.potmig.mcse

## Preparing the row names for stargazer

rownames(enter.potmig.mcse)[2] <- "Treat"
rownames(enter.potmig.mcse)[3] <- "Authorities Unfair"
rownames(enter.potmig.mcse)[4] <- "Break Law for Family"
rownames(enter.potmig.mcse)[5] <- "Risk"
rownames(enter.potmig.mcse)[6] <- "Punishment"
rownames(enter.potmig.mcse)[7] <- "Gender"
rownames(enter.potmig.mcse)[8] <- "Age"
rownames(enter.potmig.mcse)[9] <- "Education"
rownames(enter.potmig.mcse)[10] <- "Income"
rownames(enter.potmig.mcse)[11] <- "Networks Abroad"
rownames(enter.potmig.mcse)[25] <- "Treat x Authorities Unfair"
rownames(enter.potmig.mcse)[26] <- "Treat x Break Law for Family"
rownames(enter.potmig.mcse)[27] <- "Treat x Risk"
rownames(enter.potmig.mcse)[28] <- "Treat x Punishment"
rownames(enter.potmig.mcse)[29] <- "Treat x Gender"
rownames(enter.potmig.mcse)[30] <- "Treat x Age"
rownames(enter.potmig.mcse)[31] <- "Treat x Education"
rownames(enter.potmig.mcse)[32] <- "Treat x Income"
rownames(enter.potmig.mcse)[33] <- "Treat x Networks Abroad"

### Model 2: work (semi-legal migration) 

## running the model
work.potmig <- lm(worklist.pre ~ worklist.pre.dum  + fair.jam + norm.family.rev
                  + worklist.pre.dum*fair.jam +  worklist.pre.dum*norm.family.rev + risk.work.dest
                  + risk.work.dest*worklist.pre.dum + punishment.work + punishment.work*worklist.pre.dum 
                  + gender + gender*worklist.pre.dum + age + age*worklist.pre.dum 
                  + education + education*worklist.pre.dum + income.imp + income.imp*worklist.pre.dum 
                  + freq.help.net + freq.help.net*worklist.pre.dum + factor(parish), data= potmig)
summary(work.potmig)

## generate clustered standard errors
work.potmig.clust<-cluster.vcov(work.potmig, potmig$ed.parish)
work.potmig.mcse <- coeftest(work.potmig, work.potmig.clust)
work.potmig.mcse

## Preparing the row names for stargazer

rownames(work.potmig.mcse)[2] <- "Treat"
rownames(work.potmig.mcse)[3] <- "Authorities Unfair"
rownames(work.potmig.mcse)[4] <- "Break Law for Family"
rownames(work.potmig.mcse)[5] <- "Risk"
rownames(work.potmig.mcse)[6] <- "Punishment"
rownames(work.potmig.mcse)[7] <- "Gender"
rownames(work.potmig.mcse)[8] <- "Age"
rownames(work.potmig.mcse)[9] <- "Education"
rownames(work.potmig.mcse)[10] <- "Income"
rownames(work.potmig.mcse)[11] <- "Networks Abroad"
rownames(work.potmig.mcse)[25] <- "Treat x Authorities Unfair"
rownames(work.potmig.mcse)[26] <- "Treat x Break Law for Family"
rownames(work.potmig.mcse)[27] <- "Treat x Risk"
rownames(work.potmig.mcse)[28] <- "Treat x Punishment"
rownames(work.potmig.mcse)[29] <- "Treat x Gender"
rownames(work.potmig.mcse)[30] <- "Treat x Age"
rownames(work.potmig.mcse)[31] <- "Treat x Education"
rownames(work.potmig.mcse)[32] <- "Treat x Income"
rownames(work.potmig.mcse)[33] <- "Treat x Networks Abroad"

## Making the table
stargazer(enter.potmig.mcse, work.potmig.mcse, align=TRUE,
          omit=c("parish"), column.labels=c("Illegal", "Semi-Legal"), 
          omit.stat=c("LL","ser","f", "n"),
          title="Full Models, Enter",single.row=FALSE)

#################################################################
##LAW BREAKING AND LAW BENDING:##################################
## HOW INTERNATIONAL MIGRANTS NEGOTIATE WITH STATE BORDERS#######
## Schwartz, Simon, Hudson, and Johnson## #######################
## October 2020## ###############################################
### Replication Script Supplementary Information#################
#################################################################

##### PACKAGES #####
library(multiwayvcov)
library(lmtest)
library(dplyr)
library(ggplot2)

##### Data: isq #####
## Load CAT data (full)
isq <- read.csv("isq.csv", header=TRUE)


#Note: All analyses include fixed effects for parish and clustered standard errors for enumeration districts. 
#This was done to account for our sampling design, which blocked by parish and clustered by ED


### Creating the aspiring migrant subsample ###
potmig <- subset(isq, aspire.abroad.pre>4)


### Table C1: aspiring migrant models without imputed income ###

## Model 1: enter (illegal migration)
# running the model
enter.ate.noimp.potmig <- lm(enterlist.pre ~ enterlist.pre.dum + gender + age + education + income.usd + freq.help.net +  factor(parish), data=potmig)
summary(enter.ate.noimp.potmig)

# generating clustered standard errors
enter.ate.noimp.potmig.clust<-cluster.vcov(enter.ate.noimp.potmig, potmig$ed.parish)
enter.ate.noimp.potmig.mcse <- coeftest(enter.ate.noimp.potmig, enter.ate.noimp.potmig.clust)
enter.ate.noimp.potmig.mcse

## Model 2: work (semi-legal migration)
# running the model
work.ate.noimp.potmig <- lm(worklist.pre ~ worklist.pre.dum + gender + age + education + income.usd + freq.help.net +  factor(parish), data=potmig)
summary(work.ate.noimp.potmig)

# generating clustered standard errors
work.ate.noimp.potmig.clust<-cluster.vcov(work.ate.noimp.potmig, potmig$ed.parish)
work.ate.noimp.potmig.mcse <- coeftest(work.ate.noimp.potmig, work.ate.noimp.potmig.clust)
work.ate.noimp.potmig.mcse

### Figure D.1: replicating figures 1-4 from body with full sample ###

### Family and enter (illegal migration)

## running the model
enter.family <- lm(enterlist.pre ~ enterlist.pre.dum + norm.family.rev + enterlist.pre.dum*norm.family.rev 
                   + gender + age + education + income.imp + freq.help.net + factor(parish), 
                   data= isq)

summary(enter.family)

## generating clustered standard errors
enter.family.clust<-cluster.vcov(enter.family, isq$ed.parish)
enter.family.mcse <- coeftest(enter.family, enter.family.clust)

## Pull out data for marginal effects plots
# coefficients
enter.family.b2 <- enter.family.mcse[2,1]
enter.family.b3 <- enter.family.mcse[3,1]
enter.family.b22 <- enter.family.mcse[22,1] # interaction term
# variances
enter.family.varb2 <- enter.family.clust[2,2]
enter.family.varb3 <- enter.family.clust[3,3]
enter.family.varb22 <- enter.family.clust[22,22]
# covariances
enter.family.covb2b22 <- enter.family.clust[2,22]
enter.family.covb3b22 <- enter.family.clust[3,22]

list(enter.family.b2, enter.family.b3, enter.family.b22, enter.family.varb2, 
     enter.family.varb3, enter.family.varb22, enter.family.covb2b22, enter.family.covb3b22)

## transforming dataframe to graph    
enter.family.df <- isq %>% mutate(norm.family.rev = seq(1, 7, length.out=1166),
                                  enter.family.conbx = enter.family.b2+enter.family.b22*norm.family.rev,
                                  enter.family.consx = sqrt(enter.family.varb2 + enter.family.varb22*(norm.family.rev^2)+
                                                              2*enter.family.covb2b22*norm.family.rev),
                                  enter.family.ax = 1.96*enter.family.consx,
                                  enter.family.upperx = enter.family.conbx + enter.family.ax,
                                  enter.family.lowerx = enter.family.conbx - enter.family.ax)


### Family and work (semi-legal migration)

## running the model
work.family <- lm(worklist.pre ~ worklist.pre.dum + norm.family.rev + worklist.pre.dum*norm.family.rev 
                  + gender + age + education + income.imp + freq.help.net + factor(parish), 
                  data= isq)

summary(work.family)

## generating clustered standard errors
work.family.clust<-cluster.vcov(work.family, isq$ed.parish)
work.family.mcse <- coeftest(work.family, work.family.clust)

## Pull out data for marginal effects plots
# coefficients
work.family.b2 <- work.family.mcse[2,1]
work.family.b3 <- work.family.mcse[3,1]
work.family.b22 <- work.family.mcse[22,1] # interaction term
# variances
work.family.varb2 <- work.family.clust[2,2]
work.family.varb3 <- work.family.clust[3,3]
work.family.varb22 <- work.family.clust[22,22]
# covariances
work.family.covb2b22 <- work.family.clust[2,22]
work.family.covb3b22 <- work.family.clust[3,22]

list(work.family.b2, work.family.b3, work.family.b22, work.family.varb2, 
     work.family.varb3, work.family.varb22, work.family.covb2b22, work.family.covb3b22)

## transforming dataframe to graph    
work.family.df <- isq %>% mutate(norm.family.rev = seq(1, 7, length.out=1166),
                                 work.family.conbx = work.family.b2+work.family.b22*norm.family.rev,
                                 work.family.consx = sqrt(work.family.varb2 + work.family.varb22*(norm.family.rev^2)+
                                                            2*work.family.covb2b22*norm.family.rev),
                                 work.family.ax = 1.96*work.family.consx,
                                 work.family.upperx = work.family.conbx + work.family.ax,
                                 work.family.lowerx = work.family.conbx - work.family.ax)


### interaction plot with confidence intervals based on clustered SEs

# The rug
family.temp <- isq$norm.family.rev
a<-runif(1166,min=0,max=0.6)
a<-a-0.2
family.jitter<-family.temp+a

# The plot
family.df <- cbind(work.family.df, enter.family.df)
colnames(family.df) <- make.unique(names(family.df))

family.plot.combined <- ggplot(data = family.df) + 
  theme(text = element_text(size=30)) +
  scale_y_continuous("Estimated Support") + 
  scale_x_continuous(limits = c(1,7)) +
  geom_ribbon(aes(x=norm.family.rev, ymin=work.family.lowerx, ymax = work.family.upperx), fill="blue", alpha=.2) + 
  geom_line(aes(x=norm.family.rev, y=work.family.conbx, color = "Semi-Legal"), size = 1) + 
  geom_ribbon(aes(x=norm.family.rev, ymin=enter.family.lowerx, ymax = enter.family.upperx), fill="red", alpha=.2) + 
  geom_line(aes(x=norm.family.rev, y=enter.family.conbx, color = "Illegal"), size = 1)  +  
  geom_hline(yintercept = 0) +
  geom_rug(aes(x=family.jitter)) + 
  xlab("OK to Break the Law to Help Family") +
  ggtitle("Full Sample Support\n for Unauthorized Strategies") + theme(plot.title = element_text(face="bold", hjust = 0.5))+
  scale_colour_manual(name='',values=c('Semi-Legal'='blue','Illegal'='red'))
family.plot.combined

### Unfairness and enter (illegal migration)

## running the model

enter.fair <- lm(enterlist.pre ~ enterlist.pre.dum + fair.jam + enterlist.pre.dum*fair.jam 
                 + gender + age + education + income.imp + freq.help.net + factor(parish), 
                 data= isq)

summary(enter.fair)

## generating clustered standard errors
enter.fair.clust<-cluster.vcov(enter.fair, isq$ed.parish)
enter.fair.mcse <- coeftest(enter.fair, enter.fair.clust)

## Pull out data for marginal effects plots
# coefficients
enter.fair.b2 <- enter.fair.mcse[2,1]
enter.fair.b3 <- enter.fair.mcse[3,1]
enter.fair.b22 <- enter.fair.mcse[22,1] # interaction term
# variances
enter.fair.varb2 <- enter.fair.clust[2,2]
enter.fair.varb3 <- enter.fair.clust[3,3]
enter.fair.varb22 <- enter.fair.clust[22,22]
# covariances
enter.fair.covb2b22 <- enter.fair.clust[2,22]
enter.fair.covb3b22 <- enter.fair.clust[3,22]

list(enter.fair.b2, enter.fair.b3, enter.fair.b22, enter.fair.varb2, 
     enter.fair.varb3, enter.fair.varb22, enter.fair.covb2b22, enter.fair.covb3b22)

## transforming dataframe to graph    
enter.fair.df <- isq %>% mutate(fair.enter = seq(1, 7, length.out=1166),
                                enter.fair.conbx = enter.fair.b2+enter.fair.b22*fair.jam,
                                enter.fair.consx = sqrt(enter.fair.varb2 + enter.fair.varb22*(fair.jam^2)+
                                                          2*enter.fair.covb2b22*fair.jam),
                                enter.fair.ax = 1.96*enter.fair.consx,
                                enter.fair.upperx = enter.fair.conbx + enter.fair.ax,
                                enter.fair.lowerx = enter.fair.conbx - enter.fair.ax)

## Unfairness and work (semi-legal migration)
## running the model

work.fair <- lm(worklist.pre ~ worklist.pre.dum + fair.jam + worklist.pre.dum*fair.jam
                + gender + age + education + income.imp + freq.help.net + factor(parish), 
                data= isq)

summary(work.fair)

## generating clustered standard errors
work.fair.clust<-cluster.vcov(work.fair, isq$ed.parish)
work.fair.mcse <- coeftest(work.fair, work.fair.clust)

## Pull out data for marginal effects plots
# coefficients
work.fair.b2 <- work.fair.mcse[2,1]
work.fair.b3 <- work.fair.mcse[3,1]
work.fair.b22 <- work.fair.mcse[22,1] # interaction term
# variances
work.fair.varb2 <- work.fair.clust[2,2]
work.fair.varb3 <- work.fair.clust[3,3]
work.fair.varb22 <- work.fair.clust[22,22]
# covariances
work.fair.covb2b22 <- work.fair.clust[2,22]
work.fair.covb3b22 <- work.fair.clust[3,22]

list(work.fair.b2, work.fair.b3, work.fair.b22, work.fair.varb2, 
     work.fair.varb3, work.fair.varb22, work.fair.covb2b22, work.fair.covb3b22)

## transforming dataframe to graph    
work.fair.df <- isq %>% mutate(fair.work = seq(1, 7, length.out=1166),
                               work.fair.conbx = work.fair.b2+work.fair.b22*fair.jam,
                               work.fair.consx = sqrt(work.fair.varb2 + work.fair.varb22*(fair.jam^2)+
                                                        2*work.fair.covb2b22*fair.jam),
                               work.fair.ax = 1.96*work.fair.consx,
                               work.fair.upperx = work.fair.conbx + work.fair.ax,
                               work.fair.lowerx = work.fair.conbx - work.fair.ax)

### interaction plot with confidence intervals based on clustered SEs

# The rug
fair.temp <- isq$fair.jam
a<-runif(1166,min=0,max=0.6)
a<-a-0.2
fair.jitter<-fair.temp+a

# The plot
fair.df <- cbind(work.fair.df, enter.fair.df)
colnames(fair.df) <- make.unique(names(fair.df))
fair.plot.combined <- ggplot(data = fair.df) + 
  theme(text = element_text(size=30)) +
  scale_y_continuous("Estimated Support") + 
  scale_x_continuous(limits = c(1,7)) +
  geom_ribbon(aes(x=fair.jam, ymin=work.fair.lowerx, ymax = work.fair.upperx), fill="blue", alpha=.2) + 
  geom_line(aes(x=fair.jam, y=work.fair.conbx, color = "Semi-Legal"), size = 1) + 
  geom_ribbon(aes(x=fair.jam, ymin=enter.fair.lowerx, ymax = enter.fair.upperx), fill="red", alpha=.2) + 
  geom_line(aes(x=fair.jam, y=enter.fair.conbx, color = "Illegal"), size = 1)  +  
  geom_hline(yintercept = 0) +
  geom_rug(aes(x=fair.jitter),sides="b") +
  xlab("Immigration Authorities are Unfair") +
  ggtitle("Full Sample Support\n for Unauthorized Strategies") + theme(plot.title = element_text(face="bold", hjust = 0.5))+
  scale_colour_manual(name='',values=c('Semi-Legal'='blue','Illegal'='red'))

fair.plot.combined


### Perceived risk and enter (illegal migration)

## running the model
enter.risk <- lm(enterlist.pre ~ enterlist.pre.dum + risk.enter.dest + enterlist.pre.dum*risk.enter.dest 
                 + gender + age + education + income.imp + freq.help.net + factor(parish), 
                 data= isq)

summary(enter.risk)


## generating clustered standard errors
enter.risk.clust<-cluster.vcov(enter.risk, isq$ed.parish)
enter.risk.mcse <- coeftest(enter.risk, enter.risk.clust)
enter.risk.mcse

## Pull out data for marginal effects plots
# coefficients
enter.risk.b2 <- enter.risk.mcse[2,1]
enter.risk.b3 <- enter.risk.mcse[3,1]
enter.risk.b22 <- enter.risk.mcse[22,1] # interaction term
# variances
enter.risk.varb2 <- enter.risk.clust[2,2]
enter.risk.varb3 <- enter.risk.clust[3,3]
enter.risk.varb22 <- enter.risk.clust[22,22]
# covariances
enter.risk.covb2b22 <- enter.risk.clust[2,22]
enter.risk.covb3b22 <- enter.risk.clust[3,22]

list(enter.risk.b2, enter.risk.b3, enter.risk.b22, enter.risk.varb2, 
     enter.risk.varb3, enter.risk.varb22, enter.risk.covb2b22, enter.risk.covb3b22)

## transforming dataframe to graph    
enter.risk.df <- isq %>% mutate(risk.enter.dest = seq(0, 10, length.out=1166),
                                enter.risk.conbx = enter.risk.b2+enter.risk.b22*risk.enter.dest,
                                enter.risk.consx = sqrt(enter.risk.varb2 + enter.risk.varb22*(risk.enter.dest^2)+
                                                          2*enter.risk.covb2b22*risk.enter.dest),
                                enter.risk.ax = 1.96*enter.risk.consx,
                                enter.risk.upperx = enter.risk.conbx + enter.risk.ax,
                                enter.risk.lowerx = enter.risk.conbx - enter.risk.ax)

### Perceived risk and work (semi-legal migration)

## running the model
work.risk <- lm(worklist.pre ~ worklist.pre.dum + risk.work.dest + worklist.pre.dum*risk.work.dest 
                + gender + age + education + income.imp + freq.help.net + factor(parish), 
                data= isq)

summary(work.risk)

## generating clustered standard errors
work.risk.clust<-cluster.vcov(work.risk, isq$ed.parish)
work.risk.mcse <- coeftest(work.risk, work.risk.clust)

## Pull out data for marginal effects plots
# coefficients
work.risk.b2 <- work.risk.mcse[2,1]
work.risk.b3 <- work.risk.mcse[3,1]
work.risk.b22 <- work.risk.mcse[22,1] # interaction term
# variances
work.risk.varb2 <- work.risk.clust[2,2]
work.risk.varb3 <- work.risk.clust[3,3]
work.risk.varb22 <- work.risk.clust[22,22]
# covariances
work.risk.covb2b22 <- work.risk.clust[2,22]
work.risk.covb3b22 <- work.risk.clust[3,22]

list(work.risk.b2, work.risk.b3, work.risk.b22, work.risk.varb2, 
     work.risk.varb3, work.risk.varb22, work.risk.covb2b22, work.risk.covb3b22)

## transforming dataframe to graph
work.risk.df <- isq %>% mutate(risk.work.dest = seq(0, 10, length.out=1166),
                               work.risk.conbx = work.risk.b2+work.risk.b22*risk.work.dest,
                               work.risk.consx = sqrt(work.risk.varb2 + work.risk.varb22*(risk.work.dest^2)+
                                                        2*work.risk.covb2b22*risk.work.dest),
                               work.risk.ax = 1.96*work.risk.consx,
                               work.risk.upperx = work.risk.conbx + work.risk.ax,
                               work.risk.lowerx = work.risk.conbx - work.risk.ax)

### interaction plot with confidence intervals based on clustered SEs

# The rug
enter.risk.temp <- isq$risk.enter.dest
work.risk.temp <- isq$risk.work.dest
a<-runif(1166,min=0,max=0.6)
a<-a-0.2
enter.risk.jitter<-enter.risk.temp+a
work.risk.jitter<-work.risk.temp+a

# The plot
risk.df <- cbind(work.risk.df, enter.risk.df)
colnames(risk.df) <- make.unique(names(risk.df))

risk.plot.combined <- ggplot(data = risk.df) + 
  theme(text = element_text(size=30)) +
  scale_y_continuous("Estimated Support") + 
  scale_x_continuous(limits = c(0,10), labels=c(0, 2.5, 5.0, 7.5, 10)) +
  geom_ribbon(aes(x=risk.work.dest, ymin=work.risk.lowerx, ymax = work.risk.upperx), fill="blue", alpha=.2) + 
  geom_line(aes(x=risk.work.dest, y=work.risk.conbx, color = "Semi-Legal"), size = 1) + 
  geom_ribbon(aes(x=risk.enter.dest.1, ymin=enter.risk.lowerx, ymax = enter.risk.upperx), fill="red", alpha=.2) + 
  geom_line(aes(x=risk.enter.dest.1, y=enter.risk.conbx, color = "Illegal"), size = 1)  +  
  geom_hline(yintercept = 0) +
  geom_rug(aes(x=enter.risk.jitter),sides="t", color="red") +
  geom_rug(aes(x=work.risk.jitter),sides="b", color="blue") +
  xlab("Perceived Risk") +
  ggtitle("Full Sample Support\n for Unauthorized Strategies") + theme(plot.title = element_text(face="bold", hjust = 0.5))+
  scale_colour_manual(name='',values=c('Semi-Legal'='blue','Illegal'='red'))

risk.plot.combined


### Punishment and enter (illegal migration)

## running the model
enter.punishment <- lm(enterlist.pre ~ enterlist.pre.dum + punishment.enter + enterlist.pre.dum*punishment.enter 
                       + gender + age + education + income.imp + freq.help.net + factor(parish), 
                       data= isq)

summary(enter.punishment)

## generating clustered standard errors
enter.punishment.clust<-cluster.vcov(enter.punishment, isq$ed.parish)
enter.punishment.mcse <- coeftest(enter.punishment, enter.punishment.clust)

## Pull out data for marginal effects plots
# coefficients
enter.punishment.b2 <- enter.punishment.mcse[2,1]
enter.punishment.b3 <- enter.punishment.mcse[3,1]
enter.punishment.b22 <- enter.punishment.mcse[22,1] # interaction term
# variances
enter.punishment.varb2 <- enter.punishment.clust[2,2]
enter.punishment.varb3 <- enter.punishment.clust[3,3]
enter.punishment.varb22 <- enter.punishment.clust[22,22]
# covariances
enter.punishment.covb2b22 <- enter.punishment.clust[2,22]
enter.punishment.covb3b22 <- enter.punishment.clust[3,22]

list(enter.punishment.b2, enter.punishment.b3, enter.punishment.b22, enter.punishment.varb2, 
     enter.punishment.varb3, enter.punishment.varb22, enter.punishment.covb2b22, enter.punishment.covb3b22)

## transforming dataframe to graph    
enter.punishment.df <- isq %>% mutate(punishment.enter = seq(1, 4, length.out=1166),
                                      enter.punishment.conbx = enter.punishment.b2+enter.punishment.b22*punishment.enter,
                                      enter.punishment.consx = sqrt(enter.punishment.varb2 + enter.punishment.varb22*(punishment.enter^2)+
                                                                      2*enter.punishment.covb2b22*punishment.enter),
                                      enter.punishment.ax = 1.96*enter.punishment.consx,
                                      enter.punishment.upperx = enter.punishment.conbx + enter.punishment.ax,
                                      enter.punishment.lowerx = enter.punishment.conbx - enter.punishment.ax)


### Punishment and work (semi-legal migration)

## running the model
work.punishment <- lm(worklist.pre ~ worklist.pre.dum + punishment.work + worklist.pre.dum*punishment.work 
                      + gender + age + education + income.imp + freq.help.net + factor(parish), 
                      data= isq)

summary(work.punishment)

## generating clustered standard errors
work.punishment.clust<-cluster.vcov(work.punishment, isq$ed.parish)
work.punishment.mcse <- coeftest(work.punishment, work.punishment.clust)

## Pull out data for marginal effects plots
# coefficients
work.punishment.b2 <- work.punishment.mcse[2,1]
work.punishment.b3 <- work.punishment.mcse[3,1]
work.punishment.b22 <- work.punishment.mcse[22,1] # interaction term
# variances
work.punishment.varb2 <- work.punishment.clust[2,2]
work.punishment.varb3 <- work.punishment.clust[3,3]
work.punishment.varb22 <- work.punishment.clust[22,22]
# covariances
work.punishment.covb2b22 <- work.punishment.clust[2,22]
work.punishment.covb3b22 <- work.punishment.clust[3,22]

list(work.punishment.b2, work.punishment.b3, work.punishment.b22, work.punishment.varb2, 
     work.punishment.varb3, work.punishment.varb22, work.punishment.covb2b22, work.punishment.covb3b22)

## transforming dataframe to graph    
work.punishment.df <- isq %>% mutate(punishment.work = seq(1, 4, length.out=1166),
                                     work.punishment.conbx = work.punishment.b2+work.punishment.b22*punishment.work,
                                     work.punishment.consx = sqrt(work.punishment.varb2 + work.punishment.varb22*(punishment.work^2)+
                                                                    2*work.punishment.covb2b22*punishment.work),
                                     work.punishment.ax = 1.96*work.punishment.consx,
                                     work.punishment.upperx = work.punishment.conbx + work.punishment.ax,
                                     work.punishment.lowerx = work.punishment.conbx - work.punishment.ax)

### interaction plot with confidence intervals based on clustered SEs

# The rug
enter.punishment.temp <- isq$punishment.enter
work.punishment.temp <- isq$punishment.work
a<-runif(1166,min=0,max=0.6)
a<-a-0.2
enter.punishment.jitter<-enter.punishment.temp+a
work.punishment.jitter<-work.punishment.temp+a

# The plot
punishment.df <- cbind(work.punishment.df, enter.punishment.df)
colnames(punishment.df) <- make.unique(names(punishment.df))

punishment.plot.combined <- ggplot(data = punishment.df) + 
  theme(text = element_text(size=30)) +
  scale_y_continuous("Estimated Support") + 
  scale_x_continuous(limits = c(1,4)) +
  geom_ribbon(aes(x=punishment.work, ymin=work.punishment.lowerx, ymax = work.punishment.upperx), fill="blue", alpha=.2) + 
  geom_line(aes(x=punishment.work, y=work.punishment.conbx, color = "Semi-Legal"), size = 1) + 
  geom_ribbon(aes(x=punishment.work, ymin=enter.punishment.lowerx, ymax = enter.punishment.upperx), fill="red", alpha=.2) + 
  geom_line(aes(x=punishment.work, y=enter.punishment.conbx, color = "Illegal"), size = 1)  +  
  geom_hline(yintercept = 0) +
  geom_rug(aes(x=enter.punishment.jitter),sides="t", color="red") +
  geom_rug(aes(x=work.punishment.jitter),sides="b", color="blue") +
  xlab("Perceived Severity") +
  ggtitle("Full Sample Support\n for Unauthorized Strategies") + theme(plot.title = element_text(face="bold", hjust = 0.5))+
  scale_colour_manual(name='',values=c('Semi-Legal'='blue','Illegal'='red'))

punishment.plot.combined

### Table D1: Replicating Table 2 with full sample ###

### Model 1: enter (illegal migration)
## running the model
enter.full<- lm(enterlist.pre ~ enterlist.pre.dum  + fair.jam + norm.family.rev
                + enterlist.pre.dum*fair.jam +  enterlist.pre.dum*norm.family.rev + risk.enter.dest
                + risk.enter.dest*enterlist.pre.dum + punishment.enter + punishment.enter*enterlist.pre.dum 
                + gender + gender*enterlist.pre.dum + age + age*enterlist.pre.dum 
                + education + education*enterlist.pre.dum + income.imp + income.imp*enterlist.pre.dum 
                + freq.help.net + freq.help.net*enterlist.pre.dum + factor(parish), data= isq)
summary(enter.full)

## generating clustered standard errors
enter.full.clust<-cluster.vcov(enter.full, isq$ed.parish)
enter.full.mcse <- coeftest(enter.full, enter.full.clust)
enter.full.mcse

# Model 2: work (semi-legal migration)

## running the model
work.full<- lm(worklist.pre ~ worklist.pre.dum  + fair.jam + norm.family.rev
               + worklist.pre.dum*fair.jam +  worklist.pre.dum*norm.family.rev + risk.work.dest
               + risk.work.dest*worklist.pre.dum + punishment.work + punishment.work*worklist.pre.dum 
               + gender + gender*worklist.pre.dum + age + age*worklist.pre.dum 
               + education + education*worklist.pre.dum + income.imp + income.imp*worklist.pre.dum 
               + freq.help.net + freq.help.net*worklist.pre.dum + factor(parish), data= isq)
summary(work.full)

## generating clustered standard errors
work.full.clust<-cluster.vcov(work.full, isq$ed.parish)
work.full.mcse <- coeftest(work.full, work.full.clust)
work.full.mcse


### Appendix E ###

### Prep for Appendix E1 ###

risktakers <- subset(isq, risktaker>=3)
norisktakers <- subset(isq, risktaker<3)


### Perceived risk and enter among risktakers

## running the model
enter.risk.risktakers <- lm(enterlist.pre ~ enterlist.pre.dum + risk.enter.dest + enterlist.pre.dum*risk.enter.dest 
                            + gender + age + education + income.imp + freq.help.net + factor(parish), 
                            data= risktakers)

summary(enter.risk.risktakers)

## generating clustered standard errors
enter.risk.risktakers.clust<-cluster.vcov(enter.risk.risktakers, risktakers$ed.parish)
enter.risk.risktakers.mcse <- coeftest(enter.risk.risktakers, enter.risk.risktakers.clust)

## Pull out data for marginal effects plots
# coefficients
enter.risk.risktakers.b2 <- enter.risk.risktakers.mcse[2,1]
enter.risk.risktakers.b3 <- enter.risk.risktakers.mcse[3,1]
enter.risk.risktakers.b22 <- enter.risk.risktakers.mcse[22,1] # interaction term
# variances
enter.risk.risktakers.varb2 <- enter.risk.risktakers.clust[2,2]
enter.risk.risktakers.varb3 <- enter.risk.risktakers.clust[3,3]
enter.risk.risktakers.varb22 <- enter.risk.risktakers.clust[22,22]
# covariances
enter.risk.risktakers.covb2b22 <- enter.risk.risktakers.clust[2,22]
enter.risk.risktakers.covb3b22 <- enter.risk.risktakers.clust[3,22]

list(enter.risk.risktakers.b2, enter.risk.risktakers.b3, enter.risk.risktakers.b22, enter.risk.risktakers.varb2, 
     enter.risk.risktakers.varb3, enter.risk.risktakers.varb22, enter.risk.risktakers.covb2b22, enter.risk.risktakers.covb3b22)

## transforming dataframe to graph    
enter.risk.risktakers.df <- risktakers %>% mutate(risk.enter.dest = seq(0, 10, length.out=397),
                                                  enter.risk.risktakers.conbx = enter.risk.risktakers.b2+enter.risk.risktakers.b22*risk.enter.dest,
                                                  enter.risk.risktakers.consx = sqrt(enter.risk.risktakers.varb2 + enter.risk.risktakers.varb22*(risk.enter.dest^2)+
                                                                                       2*enter.risk.risktakers.covb2b22*risk.enter.dest),
                                                  enter.risk.risktakers.ax = 1.96*enter.risk.risktakers.consx,
                                                  enter.risk.risktakers.upperx = enter.risk.risktakers.conbx + enter.risk.risktakers.ax,
                                                  enter.risk.risktakers.lowerx = enter.risk.risktakers.conbx - enter.risk.risktakers.ax)



### Perceived risk and work among risktakers

## running the model
work.risk.risktakers <- lm(worklist.pre ~ worklist.pre.dum + risk.work.dest + worklist.pre.dum*risk.work.dest 
                           + gender + age + education + income.imp + freq.help.net + factor(parish), 
                           data= risktakers)

summary(work.risk.risktakers)

## generating clustered standard errors
work.risk.risktakers.clust<-cluster.vcov(work.risk.risktakers, risktakers$ed.parish)
work.risk.risktakers.mcse <- coeftest(work.risk.risktakers, work.risk.risktakers.clust)

## Pull out data for marginal effects plots
# coefficients
work.risk.risktakers.b2 <- work.risk.risktakers.mcse[2,1]
work.risk.risktakers.b3 <- work.risk.risktakers.mcse[3,1]
work.risk.risktakers.b22 <- work.risk.risktakers.mcse[22,1] # interaction term
# variances
work.risk.risktakers.varb2 <- work.risk.risktakers.clust[2,2]
work.risk.risktakers.varb3 <- work.risk.risktakers.clust[3,3]
work.risk.risktakers.varb22 <- work.risk.risktakers.clust[22,22]
# covariances
work.risk.risktakers.covb2b22 <- work.risk.risktakers.clust[2,22]
work.risk.risktakers.covb3b22 <- work.risk.risktakers.clust[3,22]

list(work.risk.risktakers.b2, work.risk.risktakers.b3, work.risk.risktakers.b22, work.risk.risktakers.varb2, 
     work.risk.risktakers.varb3, work.risk.risktakers.varb22, work.risk.risktakers.covb2b22, work.risk.risktakers.covb3b22)

## transforming dataframe to graph    
work.risk.risktakers.df <- risktakers %>% mutate(risk.work.dest = seq(0, 10, length.out=397),
                                                 work.risk.risktakers.conbx = work.risk.risktakers.b2+work.risk.risktakers.b22*risk.work.dest,
                                                 work.risk.risktakers.consx = sqrt(work.risk.risktakers.varb2 + work.risk.risktakers.varb22*(risk.work.dest^2)+
                                                                                     2*work.risk.risktakers.covb2b22*risk.work.dest),
                                                 work.risk.risktakers.ax = 1.96*work.risk.risktakers.consx,
                                                 work.risk.risktakers.upperx = work.risk.risktakers.conbx + work.risk.risktakers.ax,
                                                 work.risk.risktakers.lowerx = work.risk.risktakers.conbx - work.risk.risktakers.ax)

### interaction plot with confidence intervals based on clustered SEs

#The rug
enter.risk.risktakers.temp <- risktakers$risk.enter.dest
work.risk.risktakers.temp <- risktakers$risk.work.dest
a<-runif(397,min=0,max=0.4)
a<-a-0.2
enter.risk.risktakers.jitter<-enter.risk.risktakers.temp+a
work.risk.risktakers.jitter<-work.risk.risktakers.temp+a


# The plot
risk.risktakers.df <- cbind(work.risk.risktakers.df, enter.risk.risktakers.df)
colnames(risk.risktakers.df) <- make.unique(names(risk.risktakers.df))

risk.risktakers.plot.combined <- ggplot(data = risk.risktakers.df) + 
  theme(text = element_text(size=30)) +
  scale_y_continuous("Estimated Support") + 
  scale_x_continuous(limits = c(0,10), labels=c(0, 2.5, 5.0, 7.5, 10)) +
  geom_ribbon(aes(x=risk.work.dest, ymin=work.risk.risktakers.lowerx, ymax = work.risk.risktakers.upperx), fill="blue", alpha=.2) + 
  geom_line(aes(x=risk.work.dest, y=work.risk.risktakers.conbx, color = "Semi-Legal"), size = 1) + 
  geom_ribbon(aes(x=risk.work.dest, ymin=enter.risk.risktakers.lowerx, ymax = enter.risk.risktakers.upperx), fill="red", alpha=.2) + 
  geom_line(aes(x=risk.work.dest, y=enter.risk.risktakers.conbx, color = "Illegal"), size = 1)  +  
  geom_hline(yintercept = 0) +
  geom_rug(aes(x=enter.risk.risktakers.jitter),sides="t", color="red") +
  geom_rug(aes(x=work.risk.risktakers.jitter),sides="b", color="blue") +
  xlab("Perceived Risk") +
  ggtitle("Only Risktakers") + theme(plot.title = element_text(face="bold", hjust = 0.5))+
  scale_colour_manual(name='',values=c('Semi-Legal'='blue','Illegal'='red'))

risk.risktakers.plot.combined


### Perceived risk and enter among non-risktakers 

## running the model
enter.risk.norisktakers <- lm(enterlist.pre ~ enterlist.pre.dum + risk.enter.dest + enterlist.pre.dum*risk.enter.dest 
                              + gender + age + education + income.imp + freq.help.net + factor(parish), 
                              data= norisktakers)

summary(enter.risk.norisktakers)

## generating clustered standard errors
enter.risk.norisktakers.clust<-cluster.vcov(enter.risk.norisktakers, norisktakers$ed.parish)
enter.risk.norisktakers.mcse <- coeftest(enter.risk.norisktakers, enter.risk.norisktakers.clust)

## Pull out data for marginal effects plots
# coefficients
enter.risk.norisktakers.b2 <- enter.risk.norisktakers.mcse[2,1]
enter.risk.norisktakers.b3 <- enter.risk.norisktakers.mcse[3,1]
enter.risk.norisktakers.b22 <- enter.risk.norisktakers.mcse[22,1] # interaction term
# variances
enter.risk.norisktakers.varb2 <- enter.risk.norisktakers.clust[2,2]
enter.risk.norisktakers.varb3 <- enter.risk.norisktakers.clust[3,3]
enter.risk.norisktakers.varb22 <- enter.risk.norisktakers.clust[22,22]
# covariances
enter.risk.norisktakers.covb2b22 <- enter.risk.norisktakers.clust[2,22]
enter.risk.norisktakers.covb3b22 <- enter.risk.norisktakers.clust[3,22]

list(enter.risk.norisktakers.b2, enter.risk.norisktakers.b3, enter.risk.norisktakers.b22, enter.risk.norisktakers.varb2, 
     enter.risk.norisktakers.varb3, enter.risk.norisktakers.varb22, enter.risk.norisktakers.covb2b22, enter.risk.norisktakers.covb3b22)

## transforming dataframe to graph    
enter.risk.norisktakers.df <- norisktakers %>% mutate(risk.enter.dest = seq(0, 10, length.out=737),
                                                      enter.risk.norisktakers.conbx = enter.risk.norisktakers.b2+enter.risk.norisktakers.b22*risk.enter.dest,
                                                      enter.risk.norisktakers.consx = sqrt(enter.risk.norisktakers.varb2 + enter.risk.norisktakers.varb22*(risk.enter.dest^2)+
                                                                                             2*enter.risk.norisktakers.covb2b22*risk.enter.dest),
                                                      enter.risk.norisktakers.ax = 1.96*enter.risk.norisktakers.consx,
                                                      enter.risk.norisktakers.upperx = enter.risk.norisktakers.conbx + enter.risk.norisktakers.ax,
                                                      enter.risk.norisktakers.lowerx = enter.risk.norisktakers.conbx - enter.risk.norisktakers.ax)



### Perceived risk and work among non-risktakers

## running the model
work.risk.norisktakers <- lm(worklist.pre ~ worklist.pre.dum + risk.work.dest + worklist.pre.dum*risk.work.dest 
                             + gender + age + education + income.imp + freq.help.net + factor(parish), 
                             data= norisktakers)

summary(work.risk.norisktakers)

## generating clustered standard errors
work.risk.norisktakers.clust<-cluster.vcov(work.risk.norisktakers, norisktakers$ed.parish)
work.risk.norisktakers.mcse <- coeftest(work.risk.norisktakers, work.risk.norisktakers.clust)

## Pull out data for marginal effects plots
# coefficients
work.risk.norisktakers.b2 <- work.risk.norisktakers.mcse[2,1]
work.risk.norisktakers.b3 <- work.risk.norisktakers.mcse[3,1]
work.risk.norisktakers.b22 <- work.risk.norisktakers.mcse[22,1] # interaction term
# variances
work.risk.norisktakers.varb2 <- work.risk.norisktakers.clust[2,2]
work.risk.norisktakers.varb3 <- work.risk.norisktakers.clust[3,3]
work.risk.norisktakers.varb22 <- work.risk.norisktakers.clust[22,22]
# covariances
work.risk.norisktakers.covb2b22 <- work.risk.norisktakers.clust[2,22]
work.risk.norisktakers.covb3b22 <- work.risk.norisktakers.clust[3,22]

list(work.risk.norisktakers.b2, work.risk.norisktakers.b3, work.risk.norisktakers.b22, work.risk.norisktakers.varb2, 
     work.risk.norisktakers.varb3, work.risk.norisktakers.varb22, work.risk.norisktakers.covb2b22, work.risk.norisktakers.covb3b22)

## transforming dataframe to graph    
work.risk.norisktakers.df <- norisktakers %>% mutate(risk.work.dest = seq(0, 10, length.out=737),
                                                     work.risk.norisktakers.conbx = work.risk.norisktakers.b2+work.risk.norisktakers.b22*risk.work.dest,
                                                     work.risk.norisktakers.consx = sqrt(work.risk.norisktakers.varb2 + work.risk.norisktakers.varb22*(risk.work.dest^2)+
                                                                                           2*work.risk.norisktakers.covb2b22*risk.work.dest),
                                                     work.risk.norisktakers.ax = 1.96*work.risk.norisktakers.consx,
                                                     work.risk.norisktakers.upperx = work.risk.norisktakers.conbx + work.risk.norisktakers.ax,
                                                     work.risk.norisktakers.lowerx = work.risk.norisktakers.conbx - work.risk.norisktakers.ax)

### interaction plot with confidence intervals based on clustered SEs

#The rug
enter.risk.norisktakers.temp <- norisktakers$risk.enter.dest
work.risk.norisktakers.temp <- norisktakers$risk.work.dest
a<-runif(737,min=0,max=0.4)
a<-a-0.2
enter.risk.norisktakers.jitter<-enter.risk.norisktakers.temp+a
work.risk.norisktakers.jitter<-work.risk.norisktakers.temp+a


# The plot
risk.norisktakers.df <- cbind(work.risk.norisktakers.df, enter.risk.norisktakers.df)
colnames(risk.norisktakers.df) <- make.unique(names(risk.norisktakers.df))

risk.norisktakers.plot.combined <- ggplot(data = risk.norisktakers.df) + 
  theme(text = element_text(size=30)) +
  scale_y_continuous("Estimated Support") + 
  scale_x_continuous(limits = c(0,10), labels=c(0, 2.5, 5.0, 7.5, 10)) +
  geom_ribbon(aes(x=risk.work.dest, ymin=work.risk.norisktakers.lowerx, ymax = work.risk.norisktakers.upperx), fill="blue", alpha=.2) + 
  geom_line(aes(x=risk.work.dest, y=work.risk.norisktakers.conbx, color = "Semi-Legal"), size = 1) + 
  geom_ribbon(aes(x=risk.work.dest, ymin=enter.risk.norisktakers.lowerx, ymax = enter.risk.norisktakers.upperx), fill="red", alpha=.2) + 
  geom_line(aes(x=risk.work.dest, y=enter.risk.norisktakers.conbx, color = "Illegal"), size = 1)  +  
  geom_hline(yintercept = 0) +
  geom_rug(aes(x=enter.risk.norisktakers.jitter),sides="t", color="red") +
  geom_rug(aes(x=work.risk.norisktakers.jitter),sides="b", color="blue") +
  xlab("Perceived Risk") +
  ggtitle("No Risktakers") + theme(plot.title = element_text(face="bold", hjust = 0.5))+
  scale_colour_manual(name='',values=c('Semi-Legal'='blue','Illegal'='red'))

risk.norisktakers.plot.combined


### Prep for Figure E2 and Table E1 ### 

potmig$risk.enter.dest.d <- potmig$risk.enter.dest
potmig$risk.enter.dest.d[which(potmig$risk.enter.dest == 10)] <- NA

potmig$risk.work.dest.d <- potmig$risk.work.dest
potmig$risk.work.dest.d[which(potmig$risk.work.dest == 10)] <- NA


### Figure E2: Perceived risk with reduced risk indicator ###

### Perceived risk (reduced indicator) and enter (E2A)

## running the model
enter.risk.potmig.d <- lm(enterlist.pre ~ enterlist.pre.dum + risk.enter.dest.d + enterlist.pre.dum*risk.enter.dest.d 
                          + gender + age + education + income.imp + freq.help.net + factor(parish), 
                          data= potmig)

summary(enter.risk.potmig.d)

## generating clustered standard errors
enter.risk.potmig.clust.d<-cluster.vcov(enter.risk.potmig.d, potmig$ed.parish)
enter.risk.potmig.mcse.d <- coeftest(enter.risk.potmig.d, enter.risk.potmig.clust.d)

## Pull out data for marginal effects plots
# coefficients
enter.risk.potmig.b2 <- enter.risk.potmig.mcse.d[2,1]
enter.risk.potmig.b3 <- enter.risk.potmig.mcse.d[3,1]
enter.risk.potmig.b22 <- enter.risk.potmig.mcse.d[22,1] # interaction term
# variances
enter.risk.potmig.varb2 <- enter.risk.potmig.clust.d[2,2]
enter.risk.potmig.varb3 <- enter.risk.potmig.clust.d[3,3]
enter.risk.potmig.varb22 <- enter.risk.potmig.clust.d[22,22]
# covariances
enter.risk.potmig.covb2b22 <- enter.risk.potmig.clust.d[2,22]
enter.risk.potmig.covb3b22 <- enter.risk.potmig.clust.d[3,22]

## transforming dataframe to graph    
d.enter.risk.potmig.df <- potmig %>% mutate(risk.enter.dest.d = seq(0, 9, length.out=802),
                                            enter.risk.potmig.conbx = enter.risk.potmig.b2+enter.risk.potmig.b22*risk.enter.dest.d,
                                            enter.risk.potmig.consx = sqrt(enter.risk.potmig.varb2 + enter.risk.potmig.varb22*(risk.enter.dest.d^2)+
                                                                             2*enter.risk.potmig.covb2b22*risk.enter.dest.d),
                                            enter.risk.potmig.ax = 1.96*enter.risk.potmig.consx,
                                            enter.risk.potmig.upperx = enter.risk.potmig.conbx + enter.risk.potmig.ax,
                                            enter.risk.potmig.lowerx = enter.risk.potmig.conbx - enter.risk.potmig.ax)

### interaction plot with confidence intervals based on clustered SEs

enter.risk.potmig.d.plot <-  ggplot(data = d.enter.risk.potmig.df) + 
  theme(text = element_text(size=30)) +
  scale_y_continuous("Estimated Support", limits=c(-.5,.6)) + 
  scale_x_continuous(limits = c(0,9)) +
  geom_ribbon(aes(x=risk.enter.dest.d, ymin=enter.risk.potmig.lowerx, ymax = enter.risk.potmig.upperx), fill="red", alpha=.2) + 
  geom_line(aes(x=risk.enter.dest.d, y=enter.risk.potmig.conbx), size = 1 , color = "red")  +  
  geom_hline(yintercept = 0) + xlab("Perceived Risk") +
  ggtitle("Aspiring Migrant Support\n for Illegal Strategy") + theme(plot.title = element_text(face="bold", hjust = 0.5))+
  theme(legend.position = "none") + xlab("Perceived Risk (0-9)")

enter.risk.potmig.d.plot

### Perceived risk (reduced indicator) and work (E2B)

## running the model
work.risk.potmig.d <- lm(worklist.pre ~ worklist.pre.dum + risk.work.dest.d + worklist.pre.dum*risk.work.dest.d 
                         + gender + age + education + income.imp + freq.help.net + factor(parish), 
                         data= potmig)

summary(work.risk.potmig.d)

## generating clustered standard errors
work.risk.potmig.clust.d<-cluster.vcov(work.risk.potmig.d, potmig$ed.parish)
work.risk.potmig.mcse.d <- coeftest(work.risk.potmig.d, work.risk.potmig.clust.d)

## Pull out data for marginal effects plots
# coefficients
work.risk.potmig.b2 <- work.risk.potmig.mcse.d[2,1]
work.risk.potmig.b3 <- work.risk.potmig.mcse.d[3,1]
work.risk.potmig.b22 <- work.risk.potmig.mcse.d[22,1] # interaction term
# variances
work.risk.potmig.varb2 <- work.risk.potmig.clust.d[2,2]
work.risk.potmig.varb3 <- work.risk.potmig.clust.d[3,3]
work.risk.potmig.varb22 <- work.risk.potmig.clust.d[22,22]
# covariances
work.risk.potmig.covb2b22 <- work.risk.potmig.clust.d[2,22]
work.risk.potmig.covb3b22 <- work.risk.potmig.clust.d[3,22]

## transforming dataframe to graph    
d.work.risk.potmig.df <- potmig %>% mutate(risk.work.dest.d = seq(0, 9, length.out=802),
                                           work.risk.potmig.conbx = work.risk.potmig.b2+work.risk.potmig.b22*risk.work.dest.d,
                                           work.risk.potmig.consx = sqrt(work.risk.potmig.varb2 + work.risk.potmig.varb22*(risk.work.dest.d^2)+
                                                                           2*work.risk.potmig.covb2b22*risk.work.dest.d),
                                           work.risk.potmig.ax = 1.96*work.risk.potmig.consx,
                                           work.risk.potmig.upperx = work.risk.potmig.conbx + work.risk.potmig.ax,
                                           work.risk.potmig.lowerx = work.risk.potmig.conbx - work.risk.potmig.ax)

### interaction plot with confidence intervals based on clustered SEs

work.risk.potmig.d.plot <-  ggplot(data = d.work.risk.potmig.df) + 
  theme(text = element_text(size=30)) +
  scale_y_continuous("Estimated Support", limits=c(-.25,.6)) + 
  scale_x_continuous(limits = c(0,9), labels=c(0, 2.5, 5.0, 7.5, 9)) +
  geom_ribbon(aes(x=risk.work.dest.d, ymin=work.risk.potmig.lowerx, ymax = work.risk.potmig.upperx), fill="blue", alpha=.2) + 
  geom_line(aes(x=risk.work.dest.d, y=work.risk.potmig.conbx), size = 1 , color = "blue")  +  
  geom_hline(yintercept = 0) + xlab("Perceived Risk") +
  ggtitle("Aspiring Migrant Support\n for Semi-legal Strategy") + theme(plot.title = element_text(face="bold", hjust = 0.5))+
  theme(legend.position = "none") + xlab("Perceived Risk (0-9)")

work.risk.potmig.d.plot

### Table E1: fully interacted models with reduced risk indicator ###

### Model 1: enter (illegal migration)

## running the model
enter.potmig.d <- lm(enterlist.pre ~ enterlist.pre.dum  + fair.jam + norm.family.rev
                     + enterlist.pre.dum*fair.jam +  enterlist.pre.dum*norm.family.rev + risk.enter.dest.d
                     + risk.enter.dest.d*enterlist.pre.dum + punishment.enter + punishment.enter*enterlist.pre.dum 
                     + gender + gender*enterlist.pre.dum + age + age*enterlist.pre.dum 
                     + education + education*enterlist.pre.dum + income.imp + income.imp*enterlist.pre.dum 
                     + freq.help.net + freq.help.net*enterlist.pre.dum + factor(parish), data= potmig)
summary(enter.potmig.d)

## generating clustered standard errors
enter.potmig.d.clust<-cluster.vcov(enter.potmig.d, potmig$ed.parish)
enter.potmig.d.mcse <- coeftest(enter.potmig.d, enter.potmig.d.clust)
enter.potmig.d.mcse

### Model 2: work (semi-legal migration)

## running the model
work.potmig.d <- lm(worklist.pre ~ worklist.pre.dum  + fair.jam + norm.family.rev
                    + worklist.pre.dum*fair.jam +  worklist.pre.dum*norm.family.rev + risk.work.dest.d
                    + risk.work.dest.d*worklist.pre.dum + punishment.work + punishment.work*worklist.pre.dum 
                    + gender + gender*worklist.pre.dum + age + age*worklist.pre.dum 
                    + education + education*worklist.pre.dum + income.imp + income.imp*worklist.pre.dum 
                    + freq.help.net + freq.help.net*worklist.pre.dum + factor(parish), data= potmig)
summary(work.potmig.d)

## generating clustered standard errors
work.potmig.d.clust<-cluster.vcov(work.potmig.d, potmig$ed.parish)
work.potmig.d.mcse <- coeftest(work.potmig.d, work.potmig.d.clust)
work.potmig.d.mcse
